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Abstract. A numerical bifurcation analysis is carried out to understand the role of gravity on 
the shear-band formation in rapid granular Couette flow. At low shear rates, there is a unique 
solution with a plug near the bottom wall and a shear-layer near the top- wall; this solution mirrors 
typical shearbanding-type profiles in earth-bound shear-cell experiments. Interestingly, a stable 
plug near the top- wall is also a solution of these equations at high shear rates; there is a multitude 
of other plugged states, with the plugs being located in an ordered fashion within the Couette 
gap. The origin of such shearbanding solutions is tied to the spontaneous symmetry-breaking 
shearbanding instabilities of the gravity-free uniform shear flow, leading to both subcritical and 
supercritical pitchfork bifurcations. In the language of singularity theory, we have established that 
this bifurcation problem admits universal unfolding of pitchfork bifurcations. 

1 Introduction 

Granular materials (e.g. sand, coal, cereals, powders, etc.) are of immense importance in 
many industrial and geological processes; most agricultural and pharmaceutical products 
are in granular form. The motion of a collection of macroscopic solid particles represents 
a granular flow |7|5j as in a vibrofluidized bed. The external driving is essential to sustain 
such flows since the particle collisions are inherently dissipative, leading to a 'continual' 
energy loss. Thus, granular flows represent a driven-nonequilibrium system. In typical 'dry' 
granular flows, the effect of interstitial fluid is neglected and the interactions between grains 
are dissipative, which, in turn, leads to a wealth of interesting behaviour |7|oj . Despite their 
practical importance and non-trivial dynamics, the current understanding of granular flows 
still remains at its infancy. 

Over the last two decades, the computer simulation studies on granular shear flows have 
unveiled many interesting dynamical features (such as cluster-formation, density-waves, 
plug formation, stress-fluctuations, etc.) of such flows |8f9j . These studies have subsequently 
stimulated many works on the associated macroscopic flows described by suitable continuum 
equations. One major goal of these theoretical works has been to test the available con- 
tinuum models from the viewpoint of predicting such pattern-formation in the rapid shear 
regime. To this end, the kinetic-theory constitutive models |12|l|2j have been routinely 
used. The general concensus that appears to have emerged is that the Navier-Stokes-level 
constitutive models are able to explain many interesting dynamical states of such rapid- 
shear flows. We remark here that the non- Newtonian rheology is also important for granular 
fluids |7|3j , but theoretical works with non-Newtonian constitutive models are still lacking. 

In the present paper, we have performed linear instability and bifurcation analyses of the 
granular Couette flow under gravity using a kinetic-theory-based Newtonian hydrodynamic 
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model. Along with a synthesis of some recent work j^j, here we provide new results on 
the effect of gravity on subcritical bifurcations which firmly establish a link between our 
results and the singularity theory of bifurcations. We have shown that the origin of shear- 
band formation in such flows is tied to an instability of the associated 'gravity-free' uniform 
shear flow via the route of the universal unfolding of pitchfork bifurcations. 

2 Hydrodynamic Model: Plane Couette Flow 

The hydrodynamic model that we have employed is based on the kinetic theory of inelastic 
hard-spheres [7j. Under the assumption of instantaneous binary collisions, the kinetic theory 
provides a set of hydrodynamic equations for the density p = p p v, the hydrodynamic 
velocity u and the granular energy T, where v is the volume fraction of particles and p p is 
their intrinsic density. Note that the granular temperature is an additional (higher- order) 
field variable, defined via T =< 5u 2 /3 > where 8u is the fluctuation velocity of a particle 
over the mean motion. At the Navier-Stokes order, the hydrodynamic equations 



represent the mass, momentum and energy balances, respectively. Here Y, is the stress 
tensor, q the energy-flux, S:V« the production of fluctuation energy due to shear-work, 
D the rate of collisional dissipation of fluctuation energy and g is the acceleration due to 
gravity. Note that the collisional dissipation rate vanishes identically for perfectly elastic 
collisions (i.e. the coefficient of restitution is e = 1). As mentioned, the third equation is a 
balance for the fluctuation kinetic energy which is needed since the transport properties of 
granular fluids depend on the granular temperature. We employ the standard Newtonian 
constitutive model of an inelastic dense gas for the stress tensor X, the energy flux q and 
the collisional dissipation rate D as detailed in [TU] . 

The plane Couette flow under gravity, as shown schematically in Fig. 1, is driven by two 
oppositely-moving walls at y = —if/2 and y = H/2. In the Cartesian coordinate system, 
u(x, y, t) and v(x, y, t) denote the streamwise and transverse velocity components in the x- 
and y-directions, respectively. Using the separation between the two walls H as the length 
scale, the velocity difference between the walls U w as the velocity scale and the inverse of 
the overall shear rate H/U w = 7 _1 as the time scale, one can obtain the non-dimensional 
equations of motion in terms of dynamical variables u, u — (u,v) and T [4J. With this 
non-dimensionalization procedure, the top and bottom walls move with velocities ±1/2, 
respectively. 

For the steady (d/dt = 0), fully developed (d/dx = 0) plane Couette flow under gravity, 
the mass balance equation is identically satisfied. The x-momentum, y-momentum and 
energy balance equations reduce to 
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Fig. 1. A schematic of the plane Couette flow of monodisperse granular materials, with the 
particle size being d 



(2b) 

H~ 2 — \ k—\ + u 1 — 1 -Z> = 0, (2c) 




respectively. Here p(u,T) = fi(u, e)T, /i(z/, T) = / 2 (y,e)vT, k(u,T) = / 4 (^e)vT and 
D(z/, T) = /^(z^, e)T 3 / 2 are the nondimensional forms of the pressure, shear viscosity, pseudo- 
thermal conductivity and the collisional dissipations, respectively. The explicit forms for 
these non-dimensional functions /1 — /s are detailed elsewhere PQ. 

Only four non-dimensional control parameters are needed to describe the present prob- 
lem: the average solid fraction v avi the non-dimensional Couette gap H = H/d, the co- 
efficient of restitution e, and the Froude number Fr = U w /\fg~d (or, the non-dimensional 
weighted shear rate). Note that Fr = 77^ can be written as a ratio between two time-scales, 
where 7 = U w /H is the external shear rate and = H(dfg) 1 ^ 2 is the gravitational time 
scale. The case of zero-gravity (i.e. g = 0) corresponds to Fr^ 1 = which can also be 
recovered by considering the infinite shear-rate limit (7 — > 00). 

For all our calculations, we have adopted the no-slip and zero energy-flux boundary 
conditions. (The effects of slip-velocity and non-zero energy-flux do not introduce any new 
physics |lllll4j since the reported shear-band formation is driven by a bulk-instability.) 
With these ideal boundary conditions, the gravity-free case admits a uniform shear solution 
with constant solid fraction and granular energy: 

v{y) = const. , u(y) = y, T(y) = f 2 (p,e)/f 5 (u,e). (3) 

It is easy to verify that equations (|2afJ2c|) with Fr~ l = admit the following symmetry: 

v(y) = v(-y), u(y) = -u(-y), T(y) = T(-y). (4) 

i.e. the solid fraction and granular temperature are symmetric about the mid-plane (y = 0), 
and the velocity is antisymmetric for the uniform shear flow. However, for the plane Couette 
flow with gravity (Fr^ 1 7^ 0), this center- symmetry is no longer preserved. 
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3 Shear-banding Instability: Sub- and Supercritical Pitchfork 
Bifurcations 

For instability analysis, the uniform shear flow is perturbed by infinitesimal perturbations, 
viz., 

v(x,y,t) = u + i/(x,y,t), (5a) 
u(x,y,t) =y + u'(x,y,t), (5b) 
v(x,y,t) =v'(x,y,t), (5c) 
T(x,y,t) = f 2 /f s + T'(x,y,t). (5d) 

The dynamics of perturbations is then studied by linearizing the equations of motion around 
the base flow (v, y, 0, f%/ h) via the normal-mode analysis, 

[i/, u', v', T'\ (x, y, t) = [z>, u, v, f] (y) e ik * M . (6) 

Here k x is the wavenumber for the x-direction and the quantities with hats are complex 
amplitudes of perturbations; uj = u r + iuji is the complex frequency. Instability is associated 
with positive values of uj r . 

It is known that the uniform granular shear flow is unstable to various kinds of distur- 
bances, leading to both stationary and travelling wave patterns £Q. Here we focus only on 
a special kind of stationary instability that arises due to purely streamwise disturbances 
(k x = 0), i.e. the disturbance patterns do not vary with x\ we call it shear-banding in- 
stability since its nonlinear saturation leads to shear-banding-type patterns of alternating 
bands of dilute and dense regions in the gradient direction. For such disturbances, there 
is a minimum value of solid fraction (u av ~ 0.15) above which the flow is unstable if the 
Couette gap satisfies the following relation: 

H>nn^(iy,e), (7) 

where e ) is a complicated function of density and restitution coefficient, and n = 
1, 2, ... is the mode number (which is related to the eigenf unctions of the linearized stability 
problem) . 

It is clear that the n = 1 mode is the first to become unstable at a critical value of 
the Couette gap H = H c = iri])(y, e) (for given v av and e). Beyond H > H c , the successive 
higher-order modes (n = 2, 3, . . .) take over as the most unstable mode at H = nH c . 

Treating if as a bifurcation parameter, there is now a countably infinite number of 
pitchfork bifurcations (since the least-stable eigenvalue is real), located at H = nH c . Before 
presenting bifurcation results, let us briefly introduce the concept of pitchfork bifurcations 
by considering the simplest bifurcation problem with one order parameter 

d$ 

— = B) = b$ 3 -m = -/(-$, B), (8) 
at 

where f(4>,B) is a nonlinear scalar valued function of <P, B, called the bifurcation (control) 
parameter, and b is some real parameter. The bifurcation diagram consists of the set of 
points in the £>)-plane which satisfies /(<£, B) = 0. For pitchfork (stationary) bifurca- 
tions, the number of solutions of f(<P,B) = jumps from one to three as the bifurcation 
parameter B crosses some critical value B c . More specifically, the loss of stability of the 
base solution <£> = gives birth to two new solutions, 



<P = ±JB/b 
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which are stable. These solutions, $ = ±y£>/6, are said to be bifurcated from the base 
solution <P = 0. The bifurcation is called supercritical for b > 0, and subcritical for b < 0. 
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Fig. 2. Evidence of subcritical and supercritical bifurcations in granular Couette flow: the average 
solid fraction v av is (a) 0.16, (b) 0.3; the coefficient of restitution is e = 0.8. The stable {/unstable) 
bifurcation branch is denoted by solid (dashed) line, (c-d) Variations of solid fraction for (c) the 
first-pair (n = 1) of solution branches at H = 25 and for (d) the second-pair (n = 2) of solution 
branches at H = 50 with v av = 0.3; in each panel, the solid and dashed lines represent solutions 
corresponding to lower and upper bifurcation branches in panel (b), respectively 

For the present problem, we have chosen the order parameter to be $ = i/(l/2) = u tw , 
the value of solid fraction at the top-wall, (this can, equivalently, be tied to (v tw — v av ) or 
any other suitable quantity), and the bifurcation parameter as B = H (this can also be 
tied to (H — H c ), i.e. B = B(H, e, u av )). Figure Efa) shows a typical subcritical bifurcation 
diagram for the uniform shear case in the (H, f tlt) )-plane for a low density granular fluid 
{yav = 0.16); the coefficient of restitution is e = 0.8. The horizontal line refers to uniform 
shear for which the density is uniform across the Couette gap, i.e. v tw = v av . As we increase 
the bifurcation parameter H, the uniform shear flow loses stability to the first mode (n = 1) 
at H = H c « 64.4 and two unstable solution branches emerge from this bifurcation point. 
Each of these unstable branches eventually turns over at some lower value of H = H u m 42 
to yield a stable solution branch. For the range of Couette gaps H u < H < H c , there are 
three possible solutions, with the middle solution branch being unstable. This is an example 
of hysteretic phase transition. 
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As we increase the average density, the nature of bifurcation changes from subcritical 
to supercritical as seen in Fig. ^b) for a moderate density granular fluid {y av = 0.3). For 
these parameter combinations, the bifurcation is subcritical for v av < 0.18 and supercritical 
for v av > 0.18. From each bifurcation point, two stable solution branches emerge and such 
pairs of solutions arise due to the underlying symmetry of the uniform shear flow. Note that 
all the upper bifurcation branches asymptote to v tw — > 0.65 = v max which corresponds to 
the maximum allowable solid fraction in the system, representing the random close-packing 
limit. 

For the sake of completeness, we show the density profiles of the first two branches in Fig. 
Effe) at H = 25 and 50, respectively, in Figs. Efc-cf). The solid and dashed lines in Figs. Efc-<f) 
represent solutions belonging to the lower and upper branches of Fig. |2f&) , respectively. The 
solutions emerging from H = H c , 3H C , . . . are asymmetric, having no symmetry about the 
centerline (y = 0), since the eigenfunctions for odd modes {n = 1,3, . . .) destroy the center- 
symmetry of uniform shear. On the other hand, the solutions from H = 2H C , 4H C , . . . 
are symmetric about the centerline, since the eigenfunctions for even modes (n = 2, 4, . . .) 
satisfy the symmetry Related bifurcation results on the gravity- free uniform shear flow 
with multiple plugs are detailed elsewhere [TI]. In the following, we consider the role of 
gravity on the fate of such plugged solutions. 




Fig. 3. Effect of gravity on the bifurcation diag rani with J^av — 

0.5 and e = 0.8: (a) Fr' 1 = 0; 

(b) Fr^ 1 = 10~ 3 ; (c) Fr^ 1 = 4 x 10~ 3 . The stable (unstable) bifurcation branch is denoted by 
solid (dashed) line 
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4 Gravitational Plane Couette Flow: Gravity as an Imperfection 

Figures Of a-c) display three bifurcation diagrams at different Fr" 1 , showing the effect of 
gravity on the topology of bifurcations. The average density is set to v av = 0.5, representing 
dense flows, and other parameters are as in Fig. El For the gravity-free case (Fr" 1 = 0) 
in Fig. Efa), the topology of bifurcations is similar to that in Fig. |2^6). As we increase the 
gravitational strength to Fr" 1 = 1CU 3 (Fig. Ot), the solution branches from the even and 
odd modes peels off from each bifurcation point; this is an example of imperfect bifurcation. 
Thus, gravity acts as an imperfection in the embedded problem, destroying the bifurcation 
structure of the gravity-free case. By further increasing the gravitational strength to Fr^ 1 = 
4 x 1CT 3 (Fig. Eh), the solution branches are seen to form isolas, via merging between 
two successive bifurcation branches. Interestingly, the size of this isola also decreases with 
increasing Fr' 1 which disappears completely at Fr~ l ~ 5 x 1CT 3 . This simply implies 
that gravity does not support the associated solutions. (The sequence of birth and death 
of such isolas is detailed elsewhere j3].) The point we want to emphasize is that there is 
a critical value of Fr (~ 5 x 10 -3 for the parameter combinations of Fig. |SJ) below which 
only one solution-branch survives in the bifurcation diagram that corresponds to the lower 
branch of the first mode (n = 1) which extends upto H = oo. Note that the density profile 
corresponding to this 'unique' branch has a plug near the bottom plate as shown by the 
solid line in Fig. EJi. 




Fig. 4. Effect of average density, v av , on the profiles of (a) volume fraction, (b) granular tem- 
perature and (c) streamwise velocity at Fr^ 1 = 4 x 10 -3 , H = 25 and e = 0.8. The solutions 
correspond to the 'unique' lower bifurcation branch 
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FiguresBJ a-c) show the effect of mean density on the profiles of volume fraction, granular 
temperature and velocity at Fr~ x = 4 x 1CT 3 , H = 25 and e = 0.8. The dashed, dot-dashed 
and solid lines represent solutions at v av = 0.16, 0.3 and 0.5, respectively. We observe that 
the thickness of the plug increases with increasing u av and the plug becomes compactified. 
Within the plug, both the granular temperature and the shear-rate are relatively small. 
With increasing mean density, the plug becomes almost unsheared (i.e. the local shear rate 
being close to zero) near the bottom plate. 

4.1 Connection with Universal Unfolding of Pitchfork Bifurcation 

We have established that gravity acts as an imperfection in the present problem. In the 
parlance of singularity theory j^j, the general problem of imperfect bifurcation is addressed 
by using the theory of 'universal' unfolding which we discuss in this section. 




Fig. 5. A schematic phase-diagram for the universal unfolding of pitchfork bifurcations, showing 
all possible forms of imperfect bifurcations 



To connect our bifurcation results with universal unfoldings, let us consider the normal- 
form equation for imperfect pitchfork bifurcation jH] , 

^ = <p 3 - B<P + a + (3<P 2 = F(<P,B;a,p), (9) 
dt 

where a and (3 are the imperfection parameters which vanish identically for 'perfect' pitch- 
fork bifurcations such that 

0,0) = (10) 

From the equations of motions, we can assert that these imperfection parameters a and /3 
are functions of Fr, H and u, i.e. 

a = a(Fr~ l ,H,v) (11a) 
/3 = f3(Fr- 1 ,H,iy), (lib) 
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satisfying the following properties: 

a(0,H,v) = = /3(0,H,v). (12) 

In fact, one can think of F(<P, B; a, (3) as a perturbation of the scalar valued function f(<P, B). 
This function F($, B; a, (3) is called the universal unfolding of f(<P, B) since it encompasses 
all possible forms of bifurcation diagrams that can arise due to any form of perturbations 
on the perfect pitchfork bifurcation. 

As discussed in , the universal unfolding of the pitchfork leads to four possible bifur- 
cation diagrams in the (a, /3)-plane. In particular, two curves a = and a = /3 3 /27, denoted 
by thick solid lines in Fig. divide the (a, /5)-plane into four zones, and all parameter com- 
binations in each zone represent a canonical bifurcation diagram as shown schematically in 
this figure. Comparing our bifurcation diagrams in Fig. El with these four sets, we find that 
the bifurcation diagram for zone 1 represents the unfolding of pitchfork bifurcations in our 
case. However, by redefining our order-parameter as <P = v hw = u(— 1/2), i.e. the value of 
the solid fraction at the bottom wall, and then redrawing the bifurcation diagrams of Fig. El 
we get back the canonical set for zone 2. For this case, the smooth solution branch appears 
from the upper part of the corresponding n = 1 bifurcation branch of the gravity-free case. 




30 40 50 60 70 80 30 40 50 60 70 80 

H H 

Fig. 6. Effect of gravity on the subcritical bifurcations at v av = 0.16 and e = 0.8: (a) Fr^ 1 = 0; 
(b) Fr~ l = 2 x 10~ 4 ; (c) Fr' 1 = 4 x 10~ 4 ; (d) Fr' 1 = 4.25 x 10" 4 . The stable {unstable) 
bifurcation branch is denoted by solid (dashed) line 

We have seen that the shear-banding instabilities in granular Couette flow also lead 
to subcritical pitchfork bifurcations that occur at lower mean densities (cf. Fig. |2fr). The 



10 Meheboob Alam 



effect of gravity on such subcritical bifurcations is shown in Fig. El with parameter values 
set to v av = 0.16 and e = 0.8. Note that the smooth branch in Fig. ®b) now contains two 
saddle-nodes, joined by a 'hysteretic' branch; such 'hysteretic' jumps, however, disappear 
when the gravitational strength is strong (i.e. large Fr~ l ) as seen in Fig. [Ut- From Figs.Efc) 
andlHfc?) we observe that the size of the isola diminishes with increasing Fr _1 , eventually 
disappearing at Fr -1 ~ 4.4 x 10~ 4 for this parameter combination. Hence, under strong 
gravitational strength, the surviving attractor in the phase-space is the one that corresponds 
to a solution with a plug near the bottom wall and a dilute almost uniformly sheared layer 
near the top-wall. 

Returning back to the canonical bifurcation diagrams for universal unfolding in Fig. El 
we observe that the bifurcation diagram for zone 3 is similar to that in Fig. Again 
redefining our order-parameter as <& = v^ w = v{— 1/2) and then redrawing the bifurcation 
diagrams of Fig- El we get back the canonical set for zone 4. It is now clear that the inclusion 
of gravity naturally recovers the two 'hysteretic' bifurcation diagrams for zone 3 and zone 
4 in Fig. El Thus, all possible forms of imperfect bifurcation scenarios can be realized in 
the present bifurcation problem. We can conclude that the effect of gravity on the granular 
Couette flow truely belongs to the class of the 'universal' unfolding of pitchfork bifurcations. 

4.2 Comparison with Experiments and Discussion 

In the regime of dense flows with low shear rates, the experiments and molecular dynamics 
simulations |13ll4j suggest that the shearing mainly occurs in a thin-layer near the top wall, 
with most of the material remaining almost undeformed near the bottom wall. Recall that 
in most shear-cell experiments the Froude number is of order 10 (i.e. Fr -1 ~ 0.1). This is 
precisely what we have predicted from our analysis: at large values of Fr~ l there exists a 
unique solution with a plug near the bottom wall. Such large values of Fr~ l in earth-bound 
shear-cell experiments clearly preclude the emergence of other solutions that have plugs 
near the top wall or somewhere in between. Under micro-gravity conditions (i.e. at small 
Fr -1 ), however, it would be possible to observe such multiple-plugs, 'floating' within the 
Couette gap. Thus, we can conclude that the gravitational Couette flow admits only one 
solution, having a plug near the bottom wall and a shear-band near the top wall |2|4j . 
We have established that the shearband formation in gravitational plane Couette flow is 
tied to the shearbanding instabilities in 'uniform shear flow' via the route of the universal 
unfolding of pitchfork bifurcations. 
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